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Abstract 

A  two-phase  flow  and  multi-component  mathematical  model  with  a  complete  set  of  governing  equations  valid  in  different  components  of 
a  PEM  fuel  cell  is  developed.  The  model  couples  the  flows,  species,  electrical  potential,  and  current  density  distributions  in  the  cathode  and 
anode  fluid  channels,  gas  diffusers,  catalyst  layers  and  membrane,  respectively.  The  modeling  results  of  typical  concentration  distributions  are 
presented.  The  coupling  of  oxygen  concentration,  current  density,  overpotential  and  potential  are  shown  in  the  membrane  electrode  assembly 
(MEA).  The  model  predicted  fuel  cell  polarization  curves  for  different  cathode  pressures  compared  well  with  our  experimental  data. 
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1.  Introduction 

A  PEM  fuel  cell  sandwich  consists  of  at  least  an  anode, 
a  cathode  and  a  membrane.  Hydrogen  and  oxygen  have  to 
be  transported  to  the  anode  and  cathode  catalyst  layers  as 
reactants  by  fluid  channels  and  gas  diffuser  layers  (GDLs); 
while  the  product  water  have  to  be  transported  from  the  cat¬ 
alyst  layers  out  to  the  fluid  channels  by  GDLs.  As  results  of 
electrochemical  reactions  at  catalyst  layers,  power  and  some 
waste  heat  are  also  generated.  Two-phase  flow  and  transport 
play  essential  roles  in  the  optimum  operation  of  PEM  fuel 
cells.  First,  the  most  common  polymer  membrane  used  today 
such  as  Nation®  has  to  be  hydrated  to  ensure  high  proton  con¬ 
ductivity.  Second,  a  triple  interface  (gas,  liquid  and  solid)  in 
the  catalyst  layers  is  necessary  to  ensure  high  gas  reactant 
fluxes,  high  proton  flux  and  low  electrical  and  ionic  ohmic 
losses.  Finally,  GDLs  have  to  have  enough  open  pores  to 
transport  gas  reactants  from  fluid  channels  to  catalyst  layers 
while  effectively  removing  excessive  liquid  water  from  cata¬ 
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lyst  layers  to  fluid  channels.  If  complete  flooding  occurred  in 
the  GDLs  or  catalyst  layers,  the  maximum  gas  reactant  fluxes 
will  be  extremely  low  due  to  the  low  solubilities  of  hydrogen 
and  oxygen  in  liquid  water.  On  the  other  hand,  if  the  catalyst 
layers  and  the  membrane  are  dry,  the  ionic  resistance  will  be 
too  high  to  produce  practical  usable  currents.  Therefore,  a 
balance  of  gas-liquid  two-phase  flow  is  very  critical  for  the 
operation  of  PEM  fuel  cells.  At  present  many  of  the  transport 
phenomena  inside  a  fuel  cell  cannot  be  directly  observed  or 
measured;  thus  making  mathematical  modeling  a  critical  tool 
to  understand  such  transport  phenomena. 

The  one-dimensional  models  by  Verbrugge  and  Hill  [1,2], 
Bernardi  and  Verbrugge  [3,4],  and  Springer  et  al.  [5]  provided 
good  fundamental  bases  for  further  PEM  fuel  cell  model¬ 
ing.  Springer  et  al.  [5]  presented  empirical  relations  for  such 
parameters  as  water  diffusion  coefficient,  electro-osmotic 
drag  coefficient,  water  adsorption  isotherms,  and  membrane 
conductivities,  etc.  However,  one-dimensional  models  can¬ 
not  take  into  account  of  the  effects  of  reactants  consumption 
and  products  generation.  Fuller  and  Newman  [6],  Nguyen 
and  White  [7]  developed  two-dimensional  models.  These 
models  assume  diffusion  is  the  only  mechanism  for  oxy¬ 
gen  transport  and  did  not  consider  the  interaction  of  the 
flow  with  the  species  field  in  the  channel  and  gas  diffuser. 
Instead,  they  prescribed  concentration  boundary  conditions 
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Nomenclature 

ARm  membrane  resistance  (£2  cm2) 

Av  catalyst  surface  area  per  unit  volume 
(cm2  cm-3) 

C  mass  species  fraction 

D  diffusion  coefficient  (m2  s-1) 

Em  equivalent  weight  of  ionomer  (g  equiv- 1 ) 

F  Faraday  constant  (96,493C  mol- 1 ) 

H  Henry  constant 

i0  catalyst  exchange  current  density  (A  cm-2) 

I  current  density  (A  cm-2) 

K  absolute  permeability  (m2) 

M  molecular  weight  (kg  mol  - 1 ) 

n  number  of  electrons  or  water  transport  coef¬ 

ficient  in  the  membrane  (number  of  water 
molecules  carried  per  proton) 

N  molar  flux  (mol  cm-2  s- 1 ) 

P  pressure  (Pa) 

R  universal  gas  constant  (J  mol-1  K-1) 

T  temperature  (K) 

u  velocity  vector  or  the  v-direction  velocity  com¬ 

ponent  (ms-1) 

U  voltage  (V) 

v  y-direction  velocity  component  (ms-1) 

y  dimension  through  the  MEA  (cm) 

Zi  charge  number  of  species 

[]  activity 

Greek  letters 

a  charge  transfer  coefficient  or  net  water  transfer 

coefficient  in  the  membrane 
8  thickness  of  subscripted  layer 

£  porosity 

y  multiphase  correction  factor 

0  potential 

0  potential  difference  from  equilibrium  (overpo¬ 

tential),  V 

k  ionic  conductivity  ((£2  cm)-1) 

X  water  content  in  membrane,  mol  H2O/ 

equivalent  SO3-1 

/x  fluid  dynamic  viscosity  (N  m-2) 

v  fluid  kinematic  viscosity  (m2  s-1) 

p  density  (kg  cm-3) 

§  non-dimensional  length 

r  tortuosity 

Subscripts  and  superscripts 

0  equilibrium 

a  anode 

b  bulk  properties  or  back  pressure 

c  capillary,  cathode  or  catalyst  layer 

d  dry,  porous  media  or  electro-osmotic  drag 


diff  diffusion 

eff  effective  property  which  accounts  for  porosity, 
tortuosity  and  liquid  water 
f  ionomer  film 

i  species 

g  gas  phase 

hyd  hydraulic  permeation 

k  phase  k 

1  liquid  phase 

m  ionomer  phase 

net  net  flux 

oc  open  circuit 

solid  solid  matrix  property 

v  vapor  or  volume 

w  wet  or  water 

a  species  a 


either  at  the  gas  channel/gas  diffuser  interface,  or  at  the  cata¬ 
lyst  layer/gas  diffuser  interface.  Gurau  et  al.  [8]  was  the  first  to 
use  computation  fluid  dynamics  (CFD)  in  the  PEM  fuel  cells. 
They  presented  a  unified  approach  by  coupling  the  flow  and 
transport  governing  equations  in  the  flow  channel  and  the 
gas  diffuser,  but  only  single-phase  and  incompressible  fluid 
model  was  used. 

To  accurately  represent  the  important  transport  phenom¬ 
ena  in  PEM  fuel  cells,  it  is  well  accepted  that  a  two-phase  flow 
model  is  necessary  because  both  liquid  and  gaseous  phases 
exist  under  normal  fuel  cell  operating  conditions.  Wang  [9], 
Weber  and  Newman  [10]  provided  good  reviews  for  avail¬ 
able  two-phase  flow  models  for  PEM  fuel  cells.  Of  numerous 
two-phase  flow  models  proposed,  rigorous  multiphase  mix¬ 
ture  model  is  widely  used  because  it  can  model  the  interaction 
of  gas  and  liquid  in  the  porous  medium  with  an  acceptable 
numerical  complexity.  Wang  et  al.  [11]  proposed  the  con¬ 
cept  of  threshold  current  density  and  modeled  the  species, 
velocity  field  in  an  air  cathode  in  which  both  single  phase 
and  two-phase  flow  exist.  An  independent  parallel  work  by 
You  and  Liu  [12,13]  presented  a  two-phase  flow  model  for 
the  cathode  of  PEM  fuel  cell  and  the  model  was  used  to  dis¬ 
cuss  the  detailed  two-phase  flow  species  and  velocity  field 
as  well  as  the  effects  of  major  operating  conditions  on  the 
liquid  saturations.  Mazumder  and  Cole  [14]  also  presented  a 
multiphase  mixture  model  and  obtained  the  liquid  water  dis¬ 
tribution  in  different  components  of  PEM  fuel  cell.  Berning 
and  Djilali  [15]  presented  a  two-phase  flow  model  for  GDL 
and  flow  channel  in  both  anode  and  cathode  sides;  however, 
the  MEA  was  excluded  for  simulation.  Nam  and  Kaviany  [  1 6] 
and  Pasaogullari  and  Wang’s  [17]  works  further  elucidated 
the  effects  of  fiber  structure,  porosity,  capillary  pressures  on 
the  liquid  water  flow  in  the  GDL,  which  is  only  one  compo¬ 
nent  in  a  PEM  fuel  cell  sandwich.  In  all,  most  of  the  above 
two-phase  flow  model  efforts  except  You  and  Liu  [12,13] 
focused  on  the  GDLs,  they  either  treated  the  catalyst  layer 
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as  a  thin  film  or  just  assumed  an  operating  current  density. 
Similarly,  they  either  assumed  a  net  water  transfer  transport 
across  the  polymer  membrane  or  did  not  consider  the  water 
transport  at  all. 

There  are  several  different  mechanisms  to  decide  the  lim¬ 
iting  currents  for  PEM  fuel  cells.  The  maximum  reactant 
species  (O2  or  H2)  flux  through  GDLs  with  or  without  liq¬ 
uid  water  flooding  is  one  of  them.  As  described  in  You  and 
Liu  [18],  four  different  media  are  essentially  present  for 
electrochemical  reactions  in  the  catalyst  layers.  These  four 
media  include  pores  for  the  gas  and  for  liquid  water  transfer, 
membrane  material  for  proton  transport,  carbon  particles  for 
electron  transport  and  catalyst  particles  to  facilitate  chemical 
reactions.  The  maximum  reactant  gas  fluxes,  maximum  pro¬ 
ton  flux  across  the  catalyst  layer,  enough  amounts  of  effective 
triple  interfaces  and  the  driving  forces  for  the  electrochemi¬ 
cal  reactions  (overpotentials)  are  all  possible  factors  to  limit 
the  PEM  fuel  cell  current  density.  Still,  another  limiting  fac¬ 
tor  is  the  maximum  liquid  water  flux  through  the  membrane 
because  it  is  well  known  that  a  minimum  water  flux  is  required 
for  a  given  proton  flux  across  the  membrane.  If  the  anode  side 
is  dehydrated,  protons  will  be  difficult  to  transport  from  the 
anode  side  to  the  cathode  side,  and  thus  limit  the  PEM  fuel  cell 
current  density.  A  fuel  cell  model  will  be  incomplete  with¬ 
out  considering  all  these  components  in  the  coupled  manner, 
neither  will  it  be  complete  if  the  catalyst  layer  is  treated  as  a 
thin  film  without  considering  the  two-phase  flow  in  a  detailed 
manner. 

Therefore,  a  complete  PEM  fuel  cell  model  is  neces¬ 
sary  since  the  gas,  liquid  and  proton  transport  phenomena 
in  the  various  components  of  a  fuel  cell  are  strongly  cou¬ 
pled.  The  present  work  is  the  continuation  of  previous  works 
on  the  two-phase  flow  in  the  GDL  [12,13]  and  the  detailed 
model  on  the  catalyst  layers  [18].  In  this  work,  a  two- 
phase  multi-component  mixture  model  is  used  to  describe 
the  two-phase  flow  field  in  the  flow  channel  and  the  GDL. 
The  catalyst  layer  is  simulated  as  a  region  by  a  pseudo- 
homogeneous  model,  and  a  homogeneous  model  based  on 
the  dilute  solution  theory  is  used  for  the  membrane.  All  the 


models  are  coupled  to  study  the  whole  sandwich  of  a  PEM 
fuel  cell. 

In  this  comprehensive  model,  the  fluxes  of  hydrogen  and 
oxygen  as  well  as  the  water  production  rate  are  determined  by 
the  local  current  density.  The  net  water  transport  flux  through 
the  membrane  is  determined  by  the  electro-osmotic  drag, 
the  back  diffusion  and  hydraulic  permeation  and  is  used  as 
water  flux  boundary  conditions  in  the  anode  and  cathode.  In 
turn,  the  local  current  density  is  related  to  the  local  oxygen 
concentration  and  the  overpotential,  which  is  limited  by  the 
maximum  hydrogen  and  oxygen  flux  in  the  gas  diffuser  or 
in  the  catalyst  layer,  respectively.  In  addition,  the  effects  of 
two-phase  flow  on  the  convection  and  the  diffusion  on  the 
species  transport  are  simulated. 

2.  Mathematical  model 

2.7.  The  PEM  fuel  cell  principles 

A  simplified  diagram  of  PEM  working  principles  is  shown 
in  Fig.  1.  A  PEM  fuel  cell  consists  of  an  anode,  a  membrane 
and  a  cathode.  Each  electrode  can  be  further  divided  into 
three  regions:  the  gas  channel,  the  gas  diffuser  and  the  catalyst 
layer.  Thus,  a  PEM  fuel  cell  is  divided  into  seven  regions. 

Hydrogen  from  the  anode  side  gas  channels  transfers 
through  the  gas  diffuser  to  the  anode  catalyst  sites  where 
it  is  oxidized  under  the  action  of  catalyst  according  to  the 
following  reaction 

2H2  — >  4H+  +  4e_  (1) 

The  protons  (H+)  transport  across  the  polymer  membrane 
to  the  cathode  side.  The  electrons  are  transported  away  from 
the  anode  catalyst  layer  through  the  GDL,  bipolar  plate  and 
outer  load  to  the  cathode  catalyst  layer  where  oxygen  is 
reduced  by  the  following  reaction 

02  +  4H+  +  4e“  ->  2H20  (2) 

To  enhance  proton  transport,  the  polymer  membrane  must 
be  in  a  highly  hydrated  state.  On  one  hand,  water  is  necessary 
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Fig.  1.  Schematic  illustration  of  a  PEM  fuel  cell. 
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to  prevent  the  drying  out  of  the  membrane  to  maintain  fuel 
cell  operations.  On  the  other  hand,  too  much  water  may  cause 
flooding  in  the  pores  of  the  catalyst  layer  and  the  gas  diffuser, 
and  the  reactants  (H2,  O2)  are  not  able  to  transport  to  the 
catalyst  sites  effectively,  which  will  deteriorate  the  fuel  cell 
performance. 

From  the  above  discussions,  the  flow,  mass  transport,  cur¬ 
rent,  and  potential  in  the  different  components  in  a  PEM  fuel 
cell  are  strongly  coupled.  For  example,  the  water  flux  across 
the  membrane  depends  on  the  humidification  conditions  on 
two  sides,  the  local  current  density  and  the  membrane  prop¬ 
erties.  In  turn,  the  water  flux  influences  the  species  field  in 
the  cathode  and  anode  sides.  Furthermore,  the  oxygen  con¬ 
centration  field  affects  the  local  current  density;  and  hence 
affects  the  hydrogen,  oxygen,  and  water  mass  fraction  field 
on  both  the  anode  and  cathode  sides.  Liquid  water  also  influ¬ 
ences  the  oxygen  transport.  Therefore,  a  complete,  two-phase 
and  coupled  PEM  fuel  cell  model  is  essential  for  the  water 
and  thermal  management  as  well  as  for  the  improvement  of 
the  PEM  fuel  cell  performance. 


be  simplified,  respectively  as  [12], 

e|-(pC02)  +  V-(Ko2puC°2) 
dt 

=  V  •  (epD°2VC°2)  -  V  •  (  —  C°2L 

8  Wg  8 


e^(PCN2)  +  V-(KN2puCN2) 

Ot 

=  V  •  (epD™2  VCNz)  -  V  ■  ( —  CNlL 
8  VPgSg  8 


(7) 


(8) 


e^(pCH20)  +  V.(KH2opuCH2°) 
ot 


-V  • 


SP\  D^2°S7s\ 


-  V  • 


P 


=  V  •  (spD^2°VCU2°) 


Similar  equations  can  be  derived  for  the  anode  side. 


2.2.  Two -phase  flow  and  multi-component  model  in  the 
cathode 


A  two-phase  multi-component  mixture  model  proposed 
by  Wang  and  Chen  [19]  is  used  for  the  two-phase  flow  field 
of  humidified  air  in  the  coupled  flow  channel  and  the  gas 
diffuser.  The  governing  equations  are  presented  below,  and 
the  detailed  assumptions  and  constitutive  relations  can  be 
found  in  You  and  Liu  [12]. 

The  governing  continuity  equation  is 

dp 

£  ~  +  V  •  (pu)  =  0  (3) 

ot 

For  the  two-phase  mixture  in  the  gas  channel,  the 
Navier-Stokes  equation  is  applicable 

9(pu) 

+  V  •  (puu)  =  -  VP  +  V  •  (V/xu)  (4) 

dt 

For  the  multiphase  mixture  flow  in  the  porous  gas  diffuser, 
a  generalized  Darcy  law  is  used  [20] 

3(pu)  a 

+  V  •  (puu)  =  —VP  +  V  •  (V/xu)  -  £(eu)  (5) 

dt  K 

The  species  conservation  equation  for  the  multiphase  mix¬ 
ture  is  [18] 


s^(pCa)  +  V  •  ( yapuCa )  =  V  •  (spDVCa) 
dt 


+V- 


k 


VC“)] 


-  V  • 


Ec“j* 

k 

(6) 


For  the  specific  system  in  the  PEM  fuel  cell  cathode,  the 
mass  fraction  equations  for  oxygen,  nitrogen  and  water  can 


2.3.  Governing  equations  for  electrochemical  reaction 
in  the  catalyst  layer 


A  pseudo-homogeneous  model  [18,21,22]  is  used  to  sim¬ 
ulate  the  catalyst  layer  in  this  work.  Four  different  media  are 
usually  present  for  the  proper  function  of  catalyst  layer:  a 
diffusion  path  for  reactants  and  products  transfer,  an  ionic 
conducting  medium  for  proton  transfer,  an  electrical  con¬ 
duction  medium  for  electrons,  and  a  catalytic  surface  for 
electrochemical  reaction  to  take  place. 

We  consider  the  composite  catalyst  layer  as  a  homogenous 
medium.  From  the  oxygen  mass-current  balance  [1 8],  we  get 


(10) 


The  kinetic  expression  for  the  oxygen  reduction  rate  can  be 
described  by  Butler- Volmer  equation  if  assuming  reduction 
current  is  positive 


(11) 


where  /q"  =  ZQCf+  (Cqo/  Cq^)  and  is  reference  exchange 
current  density. 

When  water  flow  velocities  in  the  catalyst  layer  are 
neglected,  from  Ohms  law  we  have, 


d0  _  Iy_ 

d y  kq  ff 


Since 


(12) 


T)  —  E  Eq  q  —  0solid  <P  (13) 

where  q  is  negative  in  the  cathode  catalyst  layer.  The  carbon 
matrix  phase  can  be  regarded  as  equal-potential;  changes  in 
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overpotential  can  be  expressed  as 


d(— rj) 
d  y 


y 


^eff 


(14) 


According  to  Fick’s  law,  the  oxygen  flux  is  related  to  the 
oxygen  concentration  by 


,f  dCt 

AT  r^eff  ^2 

Nq2  =  ~do2-^- 


(15) 


where  Cqo  is  the  molar  concentration  of  oxygen  in  the 
ionomer,  which  is  related  to  the  molar  concentration  of  oxy¬ 
gen  in  the  gas  phase  Cqo  by  the  Henry’s  law 


0)2  =  HC^ 


(16) 


where  H  is  Henry  constant. 

Df2  in  Eq.  (15)  is  the  effective  oxygen  diffusion  coeffi¬ 
cient  given  by, 


neff  _  £l'5Do2 

uo2  ~  Z 


(17) 


Four  unknowns,  Iy,  rj ,  No2 ,  Cq2  are  described  by  Eqs. 
(10),  (11),  (14)  and  (15).  Boundary  conditions  for  the  four 
equations  are, 

At  §  =  0  (the  gas  diffuser/catalyst  layer  interface) 


Iy  =  0 


Nq2  = 


h_ 

4  F 


Co2  =  4;° 

At  §  =  1  (the  catalyst  layer/membrane  interface) 

d (-rj)  _ 
d y 


(18) 


(19) 


(20) 


K 


eff 


(21) 


2.4.  Mathematical  model  for  the  membrane 

A  generalized  expression  for  species  transport  in  the  mod¬ 
erately  dilute  solution  is  [23] 

— ZiD[FC{W  f 

Ni=  *  -AVCi- ACiVlnfl^  +  Cii;  (22) 

K1 

The  first  term  on  the  RHS  is  zero  since  water  is  not  charged, 
the  second  term  is  the  diffusion  term,  the  third  term  takes  into 
account  non-ideal  solution  behavior,  which  is  neglected  here, 
and  the  final  term  represents  the  net  transport  due  to  bulk 
flow. 

The  net  water  flux  is  the  sum  of  above  water  transport  flux 
and  the  flow  induced  by  the  electro-osmotic  drag 

net  ,  Ardiff 


K  = 

w  f 


+  +  Cwv 


ndl  ^  dCw  Km  dPw 


F 


dy 


M  dy 


(23) 


The  net  water  flux  can  be  expressed  as  net  water  transport 
coefficient  [5,24,25] 


f 

a  =  —  Nw  =  nc\  nc\\ff  n^  yd 


(24) 


where 


n&  =  2.5 


A 

22 


_  P  dCw 

^diff  —  72  w  — 


^hyd  —  77  w 


I  dy 

Km  F  d  Pw 

M  I  dy 


(25) 


(26) 


(27) 


The  number  of  water  molecules  per  sulfonate  group,  X,  is 
usually  used  to  express  the  water  content  in  the  membrane. 
It  is  related  to  the  water  molar  concentration  by  [21] 


^m,w  _  Cl&k 

w  “  cdA  +  1 


(28) 


where 


Pm,d 
~P~ 

Pm.dM' 


'in 


w 


ad  = 


Cd  = 

Pw  77m 

The  diffusive  water  transport  can  be  simplified  as 

A/diff  _  _  r)  UP 
2  V  w  —  v  ^  W 

For  one-dimensional  case 

d  a 


(29) 


(30) 


(31) 


N*F  =-Dw 


W 


W 


dym, 


(32) 


w 


where  ym5W  is  the  wet  membrane  thickness  considering  the 
change  of  membrane  dimensions  with  the  water  activity 


y- w  =  (cdA  +  l)1/3yd 


(33) 


Since 


Hfm,w 

u^w 


Hrm>w 

u^w 

dA 


dX 


0};m,w  V  dA  J  \dym  w 

solving  for  the  gradient  in  A  yields 


(34) 


dA 


n: 


diff 


w 


1 


(35) 


dym,w  7)w’W  (dCw/dA) 

After  rearranging  and  substituting  Eq.  (23)  into  Eq.  (35) 


dA 


N"et  -  (. hnd/F )  -  Cwv 


W 


dym, 


w 


72w(dCw/dA) 


(36) 


The  membrane  resistance  ARm  is  calculated  by  integrating 
over  the  membrane  thickness  by  the  following  equation 


'm 


ARm  = 


dy 


W 


r0 


k(X) 


'Sm  (cdA  +  1)1/3  dyd 
'o  k(X) 


(37) 
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2.5.  Fuel  cell  performance 
Cell  voltage  is  calculated  as 

U  =  U0C-IARm-\rjck=i  (38) 

where  \rj  c\^=l  is  the  activation  overpotential  at  the  mem¬ 
brane/catalyst  layer  interface,  ARm  can  be  predicted  from  Eq. 
(37).  The  activation  and  concentration  overpotentials  at  the 
anode  are  much  smaller  than  those  at  the  cathode,  and  thus  are 
neglected  here.  The  open  circuit  voltage  can  be  determined 
from  Nernst  equation  [26]. 

,  RT  [H2]2[02] 

Uoc  =  1.23  -  0.9  x  10"3(r  -  298)  +  —  n  - J±±J± 

4  F  [H20] 

(39) 


2.6.  Boundary  conditions  and  numerical  algorithm 


On  the  anode  side,  hydrogen  and  water  fluxes  at  the  gas 
diffuser/catalyst  layer  interface  are  related  to  local  current 
density,  respectively,  as 


Ahn  = 


N w,a  —  Li 


IF 
Mw 


F 


(40) 

(41) 


On  the  cathode  side,  oxygen  and  water  fluxes  are,  respec¬ 
tively, 


(42) 


Mw(l  +  2  a)  i 

IF 


(43) 


The  mixture  vertical  velocities  at  the  gas  diffuser/catalyst 
layer  interface  for  the  anode  side  and  the  cathode  side  are, 
respectively, 


£V 


a 


y=d\+d2 


£Vc 


y=d\+d2 


A^H2  T  ^w,a 
P 

No2  T  Afw?c 

P 


y=d\  +<i2 


(44) 


(45) 


y=d\+d2 


From  the  Darcy’s  law,  the  pressure  gradient  at  the  gas 
diffuser/catalyst  layer  interface  is 


d  P  p 
d y=~KV 


(46) 


If  natural  convection  is  neglected,  the  species  concentra¬ 
tion  gradient  can  be  obtained  from, 


d(pca ) 

=  Na~  YcPevCa 
dy 


(47) 


At  the  inlets,  the  velocities  are  given.  The  standard  exit 
boundary  and  no- slip  boundary  conditions  are  used  at  the 
channel  right  outlet  and  lower  boundary,  respectively.  For 


the  species  field,  the  inlet  species  concentrations  are  known, 
and  the  species  gradient  is  zero  at  the  lower  boundary  and  the 
right  exit  boundary. 

The  governing  equations  in  the  cathode  side,  anode  side, 
catalyst  layer  and  membrane  are  solved  simultaneously  to 
ensure  the  coupling  of  the  flow,  species,  electrical  potential 
and  current  density  distribution  in  the  cathode  and  anode  fluid 
channels,  gas  diffusers,  catalyst  layers  and  the  membrane. 

The  numerical  procedure  is  given  in  the  following.  First, 
the  overpotential  at  the  catalyst  layer/membrane  interface  is 
assumed;  the  governing  equations  in  the  catalyst  layer  are 
solved  by  the  relaxation  method  to  find  the  local  current  den¬ 
sity  and  the  oxygen  mass  flux.  From  the  local  current  density, 
the  net  water  flux  across  the  membrane  is  determined.  Then, 
the  oxygen  and  water  fluxes  are  substituted  back  into  the 
cathode  side  governing  equations  as  boundary  conditions; 
meanwhile,  the  hydrogen  and  water  fluxes  are  substituted 
back  into  the  anode  side  governing  equations.  Thus,  the  flow 
and  species  field  in  the  cathode  side  are  coupled  with  the 
flow  and  species  field  in  the  anode  side,  the  electrochemical 
reaction  in  the  catalyst  layer,  and  the  water  transport  in  the 
membrane.  Therefore,  the  flow  and  species  are  coupled  with 
the  current  and  potential  in  the  catalyst  layer  and  in  the  mem¬ 
brane.  All  the  coupled  equations  are  solved  iteratively.  Once 
the  converged  solutions  are  obtained,  the  cell  average  cur¬ 
rent  density  is  obtained  as  an  arithmetic  average  of  all  local 
current  densities  along  the  flow  direction.  By  changing  the 
input  overpotential  incrementally,  different  current  densities 
and  the  polarization  curve  is  obtained. 


3.  Results  and  discussions 

By  solving  the  flow  continuity,  momentum,  and  species 
conservation  equations  in  the  cathode  and  in  the  anode, 
respectively,  along  with  the  governing  equations  in  the 

Table  1 


The  base  case  parameters  in  the  model 


Channel  length  (cm) 

7.0 

Channel  width  (cm) 

0.1 

Gas  diffuser  thickness  (cm) 

0.03 

Catalyst  layer  thickness  (cm) 

1.29  x  1CT3 

Membrane  thickness  (cm) 

1.0  x  10“2 

Gas  diffuser  porosity 

0.4 

Catalyst  layer  porosity 

0.25 

Air  velocity  (ms-1) 

0.35 

Hydrogen  velocity  (ms-1) 

0.6 

Air  humidification  temperature  (°C) 

60 

Hydrogen  humidification  temperature  (°C) 

60 

Air  pressure  (Pa) 

3.039  x  105 

Hydrogen  pressure  (Pa) 

1.013  x  105 

Fuel  cell  temperature  (°C) 

80 

Catalyst  surface  area  (cm2  cm-3) 

1.4  x  105 

Exchange  current  density  (A  cm-2) 

4.84  x  10“8 

Cathodic  transfer  coefficient 

0.52 

Anodic  transfer  coefficient 

0.54 

Dry  membrane  density  (kgm-3) 

2.16 

Ionomer  equivalent  weight 

1100 

L.  You,  H.  Liu  /Journal  of  Power  Sources  155  (2006)  219-230 


225 


Fig.  2.  Oxygen  and  water  vapor  mass  fraction  in  the  cathode  side  for  the  base  case  with  /avg 


1.09  A  cm  2  (a)  oxygen;  (b)  water  vapor. 


catalyst  layer  and  in  the  membrane,  the  performances  of  PEM 
fuel  cell  can  be  predicted.  Table  1  lists  the  main  geometric 
and  operating  parameters  for  the  base  case. 

A  systematic  numerical  test  has  been  conducted  to  ensure 
the  simulation  results  are  grid-independent. 

3.1.  Species  concentration  in  the  cathode,  anode  and 
membrane 

Fig.  2  shows  the  cathode-side  oxygen  and  water  vapor 
mass  fraction  on  the  cathode  side  at  an  overall  fuel  cell  cur¬ 
rent  density  of  1.07  A  cm-2.  The  thickness  of  the  catalyst 
layer  in  Fig.  2a  is  enlarged  10  times  for  a  clearer  illustra¬ 
tion.  The  oxygen  fraction  decreases  along  the  flow  direction 
due  to  the  oxygen  consumption  in  the  catalyst  layer.  Consid¬ 
ering  the  enlarged  scale  for  the  catalyst  layer,  it  is  obvious 
that  the  oxygen  mass  fraction  gradient  in  the  catalyst  layer  is 
much  higher  than  those  in  the  gas  diffuser  and  gas  chan¬ 


nel.  Fig.  2b  shows  the  water  vapor  mass  fraction  in  the 
“gas  phase”  mixture.  Because  the  air  humidification  tem¬ 
perature  is  less  than  the  fuel  cell  operation  temperature,  the 
inlet  air  is  unsaturated.  The  vapor  mass  fraction  increases 
along  the  channel  due  to  the  evaporation  of  liquid  water 
from  the  catalyst  layer.  When  the  air  reaches  a  saturated 
state,  liquid  water  cannot  further  evaporate  and  remains  in 
the  liquid  state.  The  line  dividing  the  saturated  and  unsatu¬ 
rated  region  is  referred  to  “saturation  line”  hereafter.  Fig.  2b 
shows  that  the  water  vapor  fraction  increases  at  the  inlet 
of  the  channel  and  becomes  constant  after  the  saturation 
line. 

A  detailed  discussion  on  the  influences  of  fuel  cell  operat¬ 
ing  temperature,  cathode  and  anode  humidification  temper¬ 
atures  on  the  two-phase  flow  can  be  found  in  You  and  Liu 
[12]. 

Fig.  3  shows  the  corresponding  water  concentration 
contour  in  the  membrane,  where  $  =  0  is  at  the  anode 
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Fig.  3.  Water  concentration  contour  in  the  membrane  for  the  base  case  with  /avg  =  1.09  A  cm 


catalyst  layer/membrane  interface  and  §  =  1  is  at  the 
membrane/cathode  catalyst  layer  interface.  The  water 
concentration  at  the  two  boundary  sides  of  the  membrane  are 
determined  by  the  water  activities  in  the  anode  side  and  in  the 
cathode  side,  respectively  [5].  The  membrane  uptakes  water 
from  the  humidified  gases  at  two  major  surfaces.  The  water 
profile  across  the  membrane  depends  on  the  electro- osmotic 
drag  coefficient,  back  diffusion  and  hydraulic  permeation. 
Corresponding  to  the  water  activities  in  the  anode  and  cath¬ 
ode  sides,  the  water  concentration  decreases  along  the  anode 
catalyst  layer/membrane  interface  (§  =  0)  and  increases  along 
the  cathode  catalyst  layer/membrane  interface.  The  total  pro¬ 
tonic  conductance  of  the  membrane  depends  on  the  water 
concentration  profiles. 


3.2.  Coupling  of  species  concentrations  current  and 
potential  across  the  MEA 

Fig.  4  shows  the  variation  of  oxygen  concentration  across 
the  MEA  at  x  =  2.5mm  for  four  different  average  fuel  cell 
current  densities.  The  thickness  of  catalyst  layer  in  this  figure 
is  also  enlarged  10  times.  The  oxygen  mass  fraction  has  the 
highest  gradient  in  the  catalyst  layer  due  to  the  highest  oxygen 
transport  resistance  and  oxygen  consumption.  The  oxygen 
mass  fraction  gradient  in  the  gas  diffuser  is  also  higher  than 
that  in  the  fluid  channel,  owing  to  higher  transport  resistance 
and  much  lower  convection  effect  in  the  porous  gas  diffuser. 
When  the  current  density  is  low,  oxygen  flux  is  low,  and  so 
is  the  gradient  of  oxygen  concentration.  With  the  increase  of 


Fig.  4.  The  variation  of  oxygen  mass  fraction  across  the  MEA  for  the  base  case  with  four  different  current  densities. 
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Fig.  5.  The  variation  of  current  density  and  overpotential  across  the  cathode  catalyst  layer  for  the  base  case  with  four  different  current  densities. 


current  density,  more  oxygen  is  consumed  and  the  oxygen 
concentration  in  the  catalyst  layer  decreases.  With  further 
increase  of  current  density  (e.g.  I-  1.23  A  cm-2),  most  oxy¬ 
gen  is  consumed  in  a  region  of  the  catalyst  layer  close  to 
the  gas  diffuser,  which  can  be  called  the  “active  region”.  The 
remaining  region  where  no  oxygen  is  present  can  be  called  the 
“inactive  region”.  Because  of  oxygen  mass-transport  limita¬ 
tion,  the  inactive  region  does  not  contribute  to  the  production 
of  current.  At  a  given  current  density,  increasing  oxygen  con¬ 
centration  and  lowering  oxygen  reduction  overpotential  are 
favorable  to  minimize  the  inactive  region.  This  can  be  real¬ 
ized  by  decreasing  liquid  volume  fraction  and  gas  diffusion 
paths  as  well  as  increasing  gas  diffuser  porosity,  stoic  ratio 
and  air  pressure. 

Fig.  5  shows  the  variation  of  proton  current  density  and 
overpotential,  respectively,  across  the  catalyst  layer  corre¬ 
sponding  to  the  oxygen  distribution  in  Fig.  4.  The  proton 


current  in  Fig.  5a  is  zero  at  the  gas  diffuser/catalyst  layer 
interface  (§  =  0).  It  increases  across  the  catalyst  layer  due  to 
oxygen  reduction  and  reaches  the  maximum  at  the  catalyst 
layer/membrane  interface  (§  =  1).  When  the  local  current  den¬ 
sity  is  small,  the  current  is  generated  almost  homogeneously 
across  the  catalyst  layer.  With  the  increase  of  current  den¬ 
sity,  the  current  generation  rate  becomes  non-uniform.  Most 
current  is  generated  in  the  “active  region”. 

Fig.  5b  shows  the  corresponding  variation  of  overpoten¬ 
tial  across  the  catalyst  layer.  The  magnitude  of  overpotential 
corresponds  to  the  current  density.  For  the  same  oxygen  con¬ 
centration,  if  the  overpotential  is  low,  the  current  density  is 
also  low  and  vice  versa.  From  the  definition,  overpotential  is 
the  potential  difference  between  the  electrolyte  phase  and  the 
solid  phase.  The  electrolyte  phase  potential  decreases  from 
the  membrane  interface  (§  =  1)  to  the  gas  diffuser  interface 
(§  =  0)  due  to  the  Ohmic  loss  of  proton  transfer  through  the 
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Fig.  6.  The  water  profile  in  the  membrane  for  different  current  densities  for  the  base  case  with  four  different  current  densities. 


membrane  material  in  the  catalyst  layer.  Because  of  good 
electrical  conductivity,  the  potential  gradient  in  the  solid  car¬ 
bon  matrix  can  be  neglected.  Thus,  the  overpotential  has 
the  same  tendency  as  the  electrolyte  phase  potential.  The 
overpotential  has  the  highest  absolute  value  at  the  catalyst 
layer/membrane  interface  and  the  lowest  absolute  value  at 
the  catalyst  layer/gas  diffuser  interface.  The  gradient  of  over¬ 
potential  is  proportional  to  the  proton  current  if  the  membrane 
proton  conductivity  is  constant.  Therefore,  overpotential  dif¬ 
ference  for  high  current  density  is  higher  than  that  for  low 
current  density. 

Fig.  6  shows  the  water  profiles  at  x  =  2.5  mm  under  four 
different  current  densities.  Here  §  =  0  and  1  have  the  same 
designations  as  in  Fig.  3.  At  low  current  densities,  no  liquid 


water  is  present  at  the  beginning  of  cathode  side,  so  the  water 
activity  in  the  cathode  side  is  lower  than  that  in  the  anode 
side.  Therefore,  water  diffusion  has  the  same  direction  as  the 
electro-osmotic  flux.  At  higher  current  densities,  more  water 
is  produced  at  cathode  and  liquid  water  is  present.  Thus,  the 
water  activity  in  the  cathode  side  is  higher  than  that  in  the 
anode  side.  As  a  result,  the  water  diffusion  has  the  oppo¬ 
site  direction  as  the  electro-osmotic  flux.  Because  the  net 
water  flux  is  higher  at  higher  current  densities,  more  water  is 
dragged  through  the  membrane  from  the  anode  side.  There¬ 
fore,  the  water  activities  in  the  anode  side  at  higher  current 
densities  are  lower  than  those  at  lower  current  densities.  In 
contrast,  the  water  activities  in  the  cathode  side  at  higher  cur¬ 
rent  densities  are  higher  than  those  at  lower  current  densities. 


Fig.  7.  The  comparison  of  model  results  with  the  experimental  data  for  different  cathode  back  pressures  at  1,  2  and  3  bar,  respectively.  Air  side:  7hyd  =  60  °C, 
Min  =  0.35  ms-1;  hydrogen  side:  P\>  =  1  bar,  7hyd  =  60  °C,  um  =  0.6  m  s_1 ;  fuel  cell  operating  temperature  Tce\\  =  80  °C. 
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3.3.  Comparison  with  experimental  data 

Fig.  7  shows  a  group  of  polarization  curves  under  different 
cathode  pressures  from  the  model  compared  with  experimen¬ 
tal  data  [27].  The  fuel  cell  operating  temperature  is  80  °C, 
the  humidification  temperature  is  60  °C  for  both  the  cathode 
and  anode,  the  pressure  at  the  anode  is  1  atm,  the  cathode 
inlet  velocity  is  0.6  ms-1,  and  cathode  pressures  are  at  1, 
2,  and  3  atm,  respectively  [27] .  It  is  shown  that  the  model¬ 
ing  results  agree  with  the  experimental  data  reasonably  well. 
With  the  increase  of  the  cathode  pressure,  both  the  oxygen 
molar  concentration  and  exchange  current  density  increase, 
which  results  in  the  increase  of  the  current  generation  rates 
at  a  given  overpotential.  Thus,  the  fuel  cell  output  voltage 
increases. 


4.  Conclusions 

A  two-dimensional,  two-phase  flow  mathematical  model 
with  a  complete  set  of  governing  equations  for  all  the  com¬ 
ponents  of  a  PEM  fuel  cell  is  developed.  This  model  couples 
the  flow,  species,  electrical  potential  and  current  density 
distributions  in  the  two  flow  channels,  two  gas  diffusion 
layers,  two  catalyst  layers,  and  the  membrane,  respectively. 
In  the  flow  channel  and  gas  diffuser,  a  two-phase  flow 
and  multi-component  model  is  used.  In  the  cathode  cata¬ 
lyst  layer,  a  pseudo-homogeneous  model  is  used.  For  water 
transport  through  the  membrane,  the  model  includes  the 
electro-osmotic  drag,  diffusion  and  hydraulic  permeation. 
The  changes  of  membrane  proton  conductivity  and  dimen¬ 
sion  with  water  content  in  the  membrane  are  also  incorporated 
in  the  model. 

The  governing  equations  on  the  cathode  and  anode  sides, 
as  well  as  those  for  the  membrane  are  coupled  in  a  two-phase 
flow  model,  thus  they  can  provide  realistic  information  on 
various  parameter  effects.  Specifically,  water  contents  in  the 
anode  and  cathode  channels,  gas  diffusion  layers  and  catalyst 
layers  are  coupled  by  the  membrane  water  transfer;  therefore, 
realistic  water  distribution  in  the  whole  fuel  cell  sandwich 
can  be  obtained.  A  comparison  of  the  predicated  polariza¬ 
tion  curves  with  the  experimental  data  showed  reasonable 
agreement. 
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